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The Gutzwiller projection of the Schwinger-boson mean-field solution of the 2-d spin- 1/2 

. antiferromagnet in a square lattice is shown to produce the optimized, parameter-free 
(N , 

RVB ground state. We get — 0.6688 J/site and 0.311 for the energy and the staggered 

O ■ 

magnetization. The spectrum of the excited states is found to be linear and gapless near k = 



0. Our calculation suggests, upon breaking of the rotational symmetry, = 2JZ r ^Jl — 7^ 



with Z r = 1.23. 
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It is widely realized that the 2-d antiferromagnetism may play a central role in the 
copper-oxide superconductors Q. A simple yet useful model for this problem is the Heisen- 
berg model in a square lattice with nearest-neighbor couplings. Although there is no 
exact solution yet for the model, the traditional spin-wave theory, which works at zero 
temperature with broken rotational symmetry, has been known to describe quite well its 
properties^. However, away from the half-filling the long-range correlations are quickly 
destroyed. In order to understand the entire behavior of high-T c materials, a truly two- 
dimensional theory with the potential of generalizing to the doped regime is highly desired. 
The Schwinger-boson approach stands as a good candidate for this purpose. It has been 
shown to produce the qualitatively or even quantitatively correct low-temperature behavior 
of the Heisenberg model |5],f§]. Yet extension to the general t — J model is rather straight- 
forward. For example, several groups have predicted using this formalism the possible 
existence of the spiral states at finite dopings||. Despite of these apparent successes, it 
must be stressed that they are all mean-field results. The particle number constraint (hav- 
ing one boson per site) was treated on average only. The constraint which greatly limits 
the physical states of the system still presents the main technical obstacle: the including 
of many unphysical states in most of the mean-field approaches may invalidate the results. 

To go beyond the mean-field consideration, we have proposed in a series of recent work||- 
H a general scheme for carrying out a complete Gutzwiller projection on Schwinger-boson 
mean-field states of the Heisenberg and the t — J models. The new approach attempts 
to treat the constraint exactly. It has been used to study the corrections, due to the 
particle number constraint, on the mean- field results, and to find suitable variational states 
for the t — J model. But only the static properties were calculated. In this letter, we 
shall study the ground state projected from the Schwinger-boson mean-field solution of 
the Heisenberg model. In particular, we shall show the followings: 1) It is the optimized 
RVB state of the 2-d antiferromagnet within an analytic self-consistency approximation. 
2) The approximation gives — 0.3352J/bond for the ground-state energy and 0.303 for the 
staggered magnetization, in exact agreement with the modified spin-wave theory (to the 
order 1/S). Furthermore, an exact Monte Carlo evaluation yields — 0.3344 J/bond and 
0.311, much closer to the best estimated values@ —0.3346(1) J/bond and 0.31(2). This 
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confirms the optimization. 3) The excited states are explicitly constructed and analyzed. 
The excitation spectrum is found to be linear and gapless near k = 0, in agreement with 
some well-established results. The Monte Carlo calculation also yields, upon breaking of 



the rotational symmetry, = 2JZ r yl — 7^ with Z r = 1.23. The present study brings 
useful new ingredients to the RVB states in the context of high-Tc theories. 

In terms of the Schwinger-boson representation, the nearest-neighbor antiferromagnetic 
Heisenberg model may be expressed as 

1 



h=jy: 



<ij> 



where 



In the above, a unitary transformation b 



-AjjAij + s 2 ] 



E blk 



2S. 



(1) 



(2) 



bjt has been performed on one 



of the sublattices (say B). The partition function may be calculated via 
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Z((3)=Tt[V g p} = Tt\[Y[P 1 ]p, 



p = exp(—(3H) 



i=i 



where Vq (Pj) stands for the Gutzwiller projection operator for the whole lattice (the ith 
site) which enforces the constraint in (H). In |7]||, we have shown that p can be replaced 
by a differential operator Pi, 
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provided that the matrix elements of p are known: Here on = (b^, are complex numbers 
and I a* >= exp(X)o-=Ti ^IcAo-) |0 > is a usual coherent state. In the Schwinger-boson mean- 
field theory, one replaces the Hamiltonian ([I]) by a mean-field one, 



H mi = E -J E [Dijtyabj* + h.c] + A E 



(4) 



<ij>,cr 



Thus (|3]) can be used to calculate the effects of the constraint. To be definite, we shall 
restrict ourselves to the 2-d model with S = 1/2. 
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(A) The projected ground state 

We now try to project out a ground-state wave function out of the mean-field solution 
of (H). This can be done by putting (3 — » oo. The calculation of < {a]}|p m f|{o:i} > is 
straightforward. We obtain||, 

< Wll/WIK.} >= UP) exp (£[< )& LA CT + ~(<V k A + c.c.)\ j , (5) 

where 

= [(AK) sinh(/?u; k ) + cosh(^k)]- 1 , = (JD k /a*) sinh(/?u; k )^ 1) , (6) 

with u; k = ^A 2 — J 2 |.D k | 2 . Here, 6 kcr , 6 kcr and -D k are, respectively, the Fourier transforms 
of b icr , bi a and D^. Taking the uniform-bond solution = D gives -D k = zD-y^ with 
7k = [cosk x +cosk y ]/2. At /3 = oo, W 7 "^ vanishes and {frJo-} and in < {o4}|p m f > 
are fully decoupled. The ground state is simply the {b\ a } part of the Gutzwiller projection. 
This yields the wave function in the coherent-state representation, 

(7) 

{fit=0} 

where Y# is the normalization constant and [we drop, hereafter, the superscript "(2)"] 

w = 1 y ^Tkexp^k-fc-rj)] 

1 zD 
= ^E W/ kexp[^k.(r J -r,)], 77 = — - 1. (8) 

Eq. (|7|) can be easily checked by inserting p = |$g >< $g| into (H). The latter will be 
employed in the following calculations. Note that the ground state so obtained is a truly 
RVB fa/pe||, the bond strength is simply Wij. But we shall see that the |$g > possesses 
long-range antiferromagnetic correlations, in agreement with other theories. 

(B) The self-consistency approximation 

The evaluation of Y N can be reduced to a generalized loop gas problem ||,|nj. A general 
procedure was outlined in J7],|8|] . But the self-consistency approximation needs to be further 
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elucidated. The effect of Pi can be most easily handled via a 4 x 4 transfer matrix Ty 
(i, j are the lattice sites). Suppose that, at a given stage, the relevant part of the prefactor 
contains gLb^ + gu;bL + g}\bii + (They are brought down by previous differentiations. 
Prefactors not in this form will start new loops). Taking P, leads to the prefactor for site 
j via 
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(9) 



iit 
V Sji 

which defines Ty. We shall refer to the four- component vectors in (Q) as the states of the 
sites. Tracing over the matrix of a close loop yields its contribution. Yjv is then the sum of 
all possible configurations of self- avoiding loops. To illustrate this, we can decompose Y/v 
into (all jfc's below are different) 



Y 



N 



E E Yv-n-i({j fc })xTr[G Ml ---G MoJ 



(10) 



The arguments of Yv_ n _i are the sites excluded. The spin-spin correlation between different 
sublattices can be similarly calculated. Substituting b\ a — > b\ a , b ia — > d/db\ a , we have 



Si • &j ► -z E ^ 6 i-^7T^7r + 7 7 E 



One needs to modify the transfer matrices at sites i and j (they are correlated). This can 
be accomplished by multiplying simultaneously the following 4x4 matrices to the states 
of i and j 

( 



db ] db ] - 





'IT 



One then can compute the modified Y^. It turns out that only configurations with % and 
j on the same loops contribute, which are —3/4 of those in Y/v, agreeing with the rules 
proposed in ||[10 . 



We now present an analytic self-consistency approximation for Y/v and < S« ■ Sj >. 



Let us first approximate Yjv/Yj 



N-n-l 



y n . This assigns a uniform weight l/y n for 



a loop of (n + 1) bonds when both sides of ([U]) are divided by Yv. The most difficult 
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part of the problem is the self-avoiding restriction. Since the system is expected to have 
long-range correlations, it may be reasonable to ignore this restriction at the first place. 
Denote the resulting y by yo. We then recover the correct y by taking into account the 
over-counting. This allows a full analytic summation over the loops in (|T0D, which leads to 
a self-consistency equation for y 

Liy I ( i2) 

Now that an appropriate self-avoiding loop can always be defined out of a general one. To 
see this, let us start a walk at a given site, say, jo- Whenever a previous site is hit, we start 
again the walk at that site and attribute the close loop to the renormalized 1/y of the site. 
The procedure can be continued till the end of the original loop. As a result, 1/y for a 
self-avoiding loop of length (n + 1) is 



1 

y 



n 

,i=0 



1 — x V loops excluding {j = 0, • • • , i — 1} 

yo 2 



V(«+i) 3 



(13) 



where (12) has been used and the last step comes from the observation that the total 
number of sites are much larger than those excluded. Eq. (|T3|) means that one can relax 
the self- avoiding condition while replacing the correct y by yo- Another way of seeing this 
is as follows. There are infinite number of self-closing loops starting at site jo when the 
constraint is relaxed. Let a single-loop contribution be L , we can define the renormalized 

one by putting the rest into 1/y. Since L + Lq H = 1, L — 2/3 and the correct L — 1, 

we obtain 1/y = 3/(2y ). 

The spin-spin correlations can be similarly considered. We relax finding sites i and j on 
the same loop to having two independent non-self-avoiding loops starting from sites i to j 
and returning from j to i. The over-counting can be corrected again by the renormalization. 
But, comparing the present case to the above one, a new over-counting appears at site j, 
apart from those included in the renormalized 1/y. The final result should be multiplied 
by 2/3 to correct it (the overlap between the two paths are not considered in the above 
discussion). 
This gives 

1 ^ exp(— ik ■ R) 



< S A (0) • S B (R) >= 



(14) 
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where Wk|k->k+ir(i,i) = ~~ Wk, which assures no bonds between same sublattices, has be 
used. The correlations between same sublattices turn out to have the same final expression 
except for the sign factor (cf. 0). 

Let us briefly review the previous analysis || that shows the long-range antiferromag- 
netic correlations. Rewrite fll^) in the form, 

where K(x) is the complete elliptic integral of the first kind. 

Consider now the limit rj = 1, so that W'(l) is very large. The equality holds only if 
y approaches W(l) from above so that the logarithmic divergence plays a role. In fact, 
y ^ W(l) + Aexp[-n Q W'(l)n/2], with A ~ o(l) being a cutoff and n = 0.607. For large 
R, the spin-spin correlations take the form ~ exp(— R/£) with 

f = (l/4)[W'(l)/A} 1/2 exp[n 7rW'(l)/4\. 



At 7] = 1, VT'(l) = oo leads to a condensation in (14j) at k = (0,0), thus long-range 



correlations. Picking up this contribution gives the staggered magnetization M s 



M s = V-<S A (0)-S B (R) > c 



Tin 

— S 0.303. 



R^oo 2 

The ground-state energy is determined by the nearest neighbor correlations. Including 
again the condensation, we find 



J0 7T v 



-0.3352J. 



Thus, the approximation recovers the conventional spin-wave theory using the Holstein- 
Primakoff transformation (expanding the square root to the order 1/S). 

The above argument for the ordering of the system can be extended to general cases. 
Let max(|Wk|) = | Wk 1 ■ It is clear that the system posses a long-range order if and only if 
[W is equivalent to W'(l) for isotropic Wk near k ; a, S — x, y] 

\dk a dk 5 J k=kQ N^l-\Wk/y \ 
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<2' 

?/o->|W ko |+0+ 



For V^jj falling exponentially at large distances, second derivatives of \Wk\ should be ev- 
erywhere well-defined. The system is thus short-range correlated. In general, singular Wk 



will leads to non-exponential decays, such as power-law decays of W%j. It may end up with 
ordered states. For example, (||) decays as R~ 2 . This result covers the general features 
found in numerical studies |J. 

(C) The proof of the optimization 

The RVB state of (|j) is, in fact, the optimized one of this class within the analytic self- 
consistency approximation. To show this, we incorporate (|l^) into (|TJ|) with a lagrangian 
multiplier X L . Let /k = Wk/y , we are led to maximize the expression 



, where we have assumed, without loss of generality, /k = /_k = /k- The solutions of /k 
read, 



Picking up the |/k| < 1 solution leads to the one given by (|8|) with rj — Uo/Xl- It is 
straightforward to show that rj — 1 minimizes the ground-state energy as expected. 

The energy can also be rigorously evaluated without relaxing the self-avoiding restric- 
tion. This is done via the Monte Carlo method (cf. below for details). It yields, on a 
48 x 48 lattice, — 0.3344J/bond and 0.311 for the energy and the staggered magnetization 
(all digits are reliable). This strongly justifies the analytic approximation. There is no 
adjustable parameter. The energy is indeed the best presented in @. Comparing these 
values to the best estimates —0.3346(1) J/bond and 0.31(2), cf. @, our approach may have 
virtually reproduced the exact ground state. 

(D) The excited states 

Our ground state is rotationally invariant. As a result, < $g|S|$c? >= 0. Therefore the 
states Sj t0l \$G >, ot = x, y, z are the excited states of the system. More explicitly, Sj )X \$G > 
in sublattice A can be written as, 



E{(7k/k-A i )/(l-/ k 2 )} 



k 




(15) 
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where Bj x flips the "spins" of site j so that different spin components are mixed. It is easy 
to see that the three branches of states are orthogonal to each others, but not among their 
own kinds. We therefore construct the Bloch states 

11 1 V* / •! m, < <f>k.,a\H\<fik a > , , 

The calculation of the excitation spectrum reduces to the evaluations of the matrix elements 
< <f>k,a\<f>i,a >, < <frk,a\Si ■ Sj\<Pi,a > % symmetry, e M = e Ky = e M = e k (but see below). 
The computations can be most easily carried out via our matrix method: we simply need 
to modify the transfer matrices at sites i,j, k, I. The procedure for the x-branch is sketched 
below. 

Consider first the four point terms in which i,j,k,l are different. The extra matrices 
at i,j (the spin sites) have been given by (0). The extra matrices at k, I (the excitation 
sites) are simply, 

U i 



B k,x, B, x 



I 



(17) 



k,i 

There are cases where some of the sites coincide. We shall classify them as three-point and 
two-point terms. The corresponding matrices can be easily deduced from ( |TT| ) and ( p]) by 
multiplications of two or more matrices at one site. One needs to pay some attentions to 
the orders of the matrices: The relevant operator to be averaged is -B^Sj • SjBj x , with 
p = \$>g >< &g\ (which provides the "loop gas"). As an example, the three-point term 
with k — i have the matrices of i,j which are simply ( |TT]) left-multiplied by the matrix of 
B k , x in (0), 

, / 





One can likewise consider the remaining cases. Inserting these matrices at i,j,k,l, one 
then computes the modified Y^. The non-zero loops involving the four sites are shown in 
Figure 1, where the values attached are the ratios of their contributions to those without 
the inserting matrices. They are used in the Monte Carlo evaluations. 

In our calculation, the loop configurations are updated by randomly choosing a pair 
of next nearest neighbor sites and exchanging their loop connections with a probability 



satisfying the detailed balance conditional . The ground-state energy per bond is obtained 
by sampling over nearest neighbor bonds on same loops (with a weight —3/4). To find 
the excited energies, we randomly pick up a spin bond i,j and the first excited site k. 
Summing the sites I over the whole lattice for a given k can be managed as follows: Find 
the loop containing k. Only the sites on this loop contribute to the denominator in fll6|) . 
This algorithm ensures the positivity of the denominator except for k = (tt, tt) where it 
is rigorously zero (thus the state is not well defined). To obtain the numerator, classify 
the location of the bond into three cases: a) Both i,j on the loop; b) Both i,j not on the 
loop; c) One of say i, on the loop. In cases a) and b), again only the sites on the loop 
of k contribute, while in case c) the sites on the loop containing j are considered. In all 
the cases, the weights in Figure 1 are used as mentioned. Figure 2 presents the numerical 
result (plotted in o) on a 10 x 10 lattice. Two points are immediately clear: (1) The gap 
at k = (0, 0) is virtually zero. (2) The spectrum is linear at small k, in agreement with 
some of the well-established results. The technical difficulty here is that one needs at least 
the precise value of the finite-size ground-state energy within a relative error of 10~ 4 . Note 
that, in contrast, the dimer state with only the nearest neighbor RVB bond always gives a 
negative energy gap at k = (beyond the error range). Thus the latter is, as a matter of 
fact, not a suitable trial ground state. 

The above spectrum does not agree in details with those obtained via more sophisticated 
spin-wave or related theories E3]. We now discuss the breaking of the rotational symmetry 
which is, we suspect, the source of the discrepancies. At zero temperature, it is generally 
believed that the system should suffer a symmetry breaking such that the spins like to 
align up on one sublattice and down on the other. When this happens, different branches 
of the excited states proposed above are no longer orthogonal. One is therefore forced to 
reconstruct the excited states. The details might be quite complicated and we shall explore 
it elsewhere. How is this going to show up in the above calculation? A natural guess is 
that, since the periodicity of the system is reduced, we only sample over terms with k, I 
belonging to same sublattices. 

We stress that this is based on our intuitive understanding of the system and has not 
been proved. The result (plotted in •) is also shown in Figure 2. It turns out to fit the 
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renormalized spin-wave result ek = 2JZ r yl — 7^, with Z r = 1.23, in excellent agreement 
with other results (cf., in particular, [IT] done on a supercomputer!). This greatly narrows 



the differences between the RVB and the spin-wave approaches, although details remain to 
be clarified. 

In conclusion, we have obtained in this work an optimized, parameter-free RVB state 
for the 2-d antiferromagnetic Heisenberg model. We have also developed effective methods 
for calculations of various physical quantities, in particular, the excitation spectrum of the 
ground state (which has been a long-standing difficulty in this context). Our results agree 
with, and in some aspects, are better than the conventional spin-wave theory. 

This work was supported by the National Science Council and the National Educational 
Committee of China. 
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Figure captions 

Fig. 1 A list of the non-zero loops involving the four sites k, I and their contributions. 

Fig. 2 The excitation spectrum e^/J, plotted in o, along the (1,0) axis (on the right-hand 
side) and the (1,1) direction. The gap at k = (0,0) is virtually zero. Plotted in • 
is 6k/J (x3) upon breaking of the rotational symmetry. Solid curve (x3), with an 
upward shift = 0.24 (suppressed as L increases), is the renormalized spin-wave result 
with Z r = 1.23. 
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Figure 1 




Figure 2 



